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Abstract. The propagation and roughening of a fluid-gas interface through a disordered medium in the case 
of capillary driven spontaneous imbibition is considered. The system is described by a conserved (model 
B) phase-field model, with the structure of the disordered medium appearing as a quenched random field 
a(x). The flow of liquid into the medium is obtained by imposing a non-equilibrium boundary condition 
on the chemical potential, which reproduces Washburn's equation H ~ t^^^ for the slowing down motion of 
the average interface position H. The interface is found to be superrough, with global roughness exponent 
X ~ 1.25, indicating anomalous scaling. The spatial extent of the roughness is determined by a length scale 
~ H^^^ arising from the conservation law. The interface advances by avalanche motion, which causes 
temporal multiscaling and qualitatively reproduces the experimental results of Horvath and Stanley [Phys. 
Rev. E 52 5166 (1995)] on the temporal scaling of the interface. 

PACS. 47.55. Mh Flows through porous media - 05.40-|-j Fluctuation phenomena, random processes, and 
Brownian motion - 68.35.Ct Interface structure and roughness 



1 Introduction 

The dynamics of driven interfaces in disordered media is 
a subject of intense interest in nonequilibrium statistical 
mechanics. It is well established that for sufhciently strong 
driving, the interface feels an effective smeared out "ther- 
mal" noise and its fluctuations present all the typical phe- 
nomena of scale invariance of driven systems jl]. In the 
opposite case of weak driving, the quenched nature of the 
noise becomes apparent and the interface may reach a 
pinned state, characterised by completely different scaling 
exponents 

An apparently easy experiment to perform is to mon- 
itor the motion of an invading liquid front in a porous 
medium. Many experiments have been done with Hele- 
Shaw cells [§,§,Hj6|^, and the spontaneous imbibition of 
water in paper has also been considered ||^,|9|, |l0|jTllJl^ , |l3| ,p^ 1 
In this last case, capillary forces arising from the porous 
structure drive the liquid until loss of water by evapo- 
ration or hydrostatic pressure balance the driving. The 
eventual pinning of the front has received a lot of at- 
tention [p|9|, [l^ , |l^ , and some experiments also examined 
the complete dynamical process [f0|O,ll^,|l^ . It is gen- 
erally believed that phenomena in the critical region of 
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the depinning transition can be described by local the- 
ories, i.e., the physics is governed by an equation that 
couples the interface locally with itself and the quenched 
randomness 1^ . The spatial configurations of a pinned im- 
bibition front seem to exhibit scaling properties well de- 
scribed by the "Directed Percolation Depinning" (DPD, or 
Quenched Kardar-Parisi-Zhang, QKPZ) universality class 
[p||9| Jl^ , |l8| . An intuitively motivated lattice model of DPD 
compares well to experimental findings on the stopped 
front 1^,^, and a modified version of the model addresses 
the influence of e vap oration on the scaling properties of 
the pinned front |12[ . 

However, these models neglect the fact that liquid has 
to be transported through the medium in order to drive 
the front, a nontrivial phenomenon in itself p9| , p0pl| ]. For 
example, viscous fluid transport explains why the invading 
Continuously slows down even without evaporation or 
a result that has been well established in the lit- 
erature [^ , p3|j2^ . In local models this has to be put in 
rather artificially jl^. The temporal correlations of the 
fluctuations should also reflect this nonlocality. It is thus 
doubtful whether any local model can explain the exper- 
imental results of Horvath and Stanley, focusing on the 
dynamical scaling of the interface . 

The main concern in this paper is to analyse in de- 
tail a simple model of a propagating liquid-gas interface 
in a disordered medium, already introduced in Ref. [p5| . 
The previous companion paper ||2^ presents a general 
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overview of the experimental and macroscopic aspects of 
imbibition and discusses in particular the role of the fluid 
conservation law. The main goal in building the present 
model is thus to incorporate local liquid conservation to 
the interfacial dynamics. In Section ^ it is shown that 
this can be achieved through a generalised Cahn-Hilliard 
equation. The model can then be applied to two differ- 
ent physical situations: a freely rising front, and a front 
propagating against a steady motion of the paper towards 
the liquid reservoir, leading to stationary fluctuations. A 
central feature of the model is the existence of a length 
scale , emerging from the interplay between interfacial 
tension and liquid conservation. This reflects the inher- 
ently non-local nature of the dynamics. The dynamical 
evolution of the model is numerically integrated in Sec- 
tion ||, for the two different setups. In Section ^ further 
implications of our modeling are discussed, and dynamical 
scaling in the experiments of Horvath and Stanley is 
reinterpreted. We conclude with suggestions for future ex- 
perimental work. The Appendices contain generalisations 
of the model to cases where gravity and/or evaporation 
might be important and a brief description of the steps 
required to obtain the interface equation, Eq. (0). 



2 Phase Field Model of Imbibition 



resulting equation of motion, 

9t<^(x,t) = VV(x,i) (2) 
= V2 [-0(x, t) + 03(x, t) - V^<f>{^, t) - a(x)] , 

is essentially the Cahn-Hilliard equation |2^, also used 
to study critical fluctuations and phase ordering in pres- 
ence of a conservation law (model B dynamics p^ , p9| ). 
The variable q;(x) here plays the role of the local chemical 
potential at the interface thus controlling the flux. 

In model B dynamics pof , the domain walls are driven 
by the difference between incoming and outgoing current 
j = — V/i. In the sharp interface limit, and for a slowly mov- 
ing front, /Lt(x, t) changes quasistatically, always satisfying 
V^/i = in the bulk, plus the appropriate boundary condi- 
tions. At the interface, /i must obey the Gibbs- Thomson 
boundary condition 

Acp fiU ^AV- aJC, (3) 

where JC is the curvature, a = 2\/2/3 is the surface ten- 
sion of the model, the miscibility gap A(j) — 0+ — 
and AV — V{4)+) — V{4)-). The quantities 0± are the 
equilibrium values of the phase fleld, defined by the usual 
tangent construction p9| , p0t . The interface motion is then 
determined by the normal velocity w„ = —d„p\^- 



2.1 Definition of the model 

A full model of the dynamics of a liquid invading a random 
medium based on a microscopic treatment is a formidable 
problem. The dynamics of the advancing interface be- 
tween the liquid and the (dry) solid should however be 
amenable to a discussion at the coarse grained level. In 
this spirit, a phase-field model is used to represent the 
spatial configuration of wet and dry "phases" . The field of 
interest is a locally conserved quantity (/)(x, t) , defined on 
the half-plane {x= {x, y)\y > 0}, with values (p = -|-1(— 1) 
for the wet (dry) phase. A free energy of the form = 
J df-x [(V0)^/2 -I- ¥{(())] is chosen, with a double weU po- 
tential V((j)) of the standard Ginzburg-Landau form, to 
which is added a linear tilt, 

l/(x,0(x,t)) = -i <^2(x,t) + i 04(x,t)_a(x)0(x,O. (1) 

The double well potential, together with the gradient term, 
ensures the existence of a well defined interface, and the 
quenched random field a(x) represents the random (coarse 
grained) structure of the medium. The first and the sec- 
ond moments associated with the underlying distribution 
of the random medium are given by (q;(x)) = a, and 
(Q(x)a(x'))— — (Z\a)^(5(x— x'). It is thus assumed that a 
is spatially uncorrelated, which may be a good approxima- 
tion in the case of ordinary paper (the areal mass density 
has only short-range correlations |M). 

The dynamics of the conserved variable 0(x, t) is de- 
termined by a continuity equation dt4> + V-j = 0, where 
the current j(x, t) = — V^(x, t) is related to the gradi- 
ent of the chemical potential /i(x, i) = ~5!F/5(j){x,t). The 



2.2 Freely rising and stationary fronts 

With appropriate boundary conditions, the model pre- 
sented above can encompass the typical experimental se- 
tups of imbibition, as depicted in Fig. |l|. These are such 
that the value of the chemical potential fi{x,y = 0) = 
ao ^ a is imposed at the bottom end while the top end 
of the system is kept dry (i.e., dyfx{y oo) = and 
4>{y —> oo) = —1). This concept can easily be explained 
by a simplified situation where the quenched field a(x) is 
homogeneous and equal to a constant a. In this case, an 
equilibrium interface would be obtained by letting /i = —a 
throughout the whole system, with (j>± ~ ±1. On the other 
hand, imposing the boundary condition /Lt(a;, y = 0) = ao 
creates an imbalance in the chemical potential causing the 
interface to advance. If the interface at time t is at a height 
H{t), the chemical potential, as given by the Laplace equa- 
tion, is 

- y 

H{x,t) = fi{y,t) = ao - (a " "o)-^^, for y < H; 

fi{y,t) = -a, for y>H, (4) 

resulting in a time evolution Q 

dH{t) a — ckq , . 

dt ~ 2H{t) ■ ^ ' 

Thus, the further the interface is from the reservoir, the 
smaller its velocity. This classical result, known as the 

^ This is valid only if a — qq 1. A more complete expression 
is given in Appendix 
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Washburn equation, is well established experimentally, al- 
though discrepancies may arise |22||2^ , p^ , pl| . 

A slight modification to the model can be used to re- 
produce the experimental setup of Horvath and Stanley 
p!3| . In their experiment, the interface is forced to stay 
at a fixed mean height H by constantly pulling down the 
paper towards the reservoir of liquid. Within the phase 
field model it is easy to add a constant downward drift 
V — — wy, so that 

^(^(x,*) + vV?!)(x,i) = 
[-0(x, t) + 03(x, t) - V20(x, t) - a(x- vt)] , (6) 

keeping the same boundary conditions as for the freely 
rising column. Thus the interface between the wet and the 
dry region is kept at a height H where a rising interface 
would have a velocity — v, or 



H=- 



"0 



2v 



2.3 Equation of motion and correlation length 



(7) 



In presence of quenched disorder, via the field q;(x), the 
interface will start to roughen, as shown in Fig. |^. An im- 
mediately noticeable feature is that the interface looks ex- 
tremely rough locally but appears smooth on large length 
scales. This should indeed be expected intuitively, since 
the physics of the phenomenon is such that parts of the 
interface ahead of the average interface position have a 
smaller instantaneous local velocity. They are thus even- 
tually "caught up" by the average interface. Likewise, re- 
tarded parts of the interface tend to catch up with the 
average interface position. This idea is indeed confirmed 
by the numerical results, and can further be used to define 
the spatial range over which the correlated roughness may 
be seen. 

A key step in understanding the physics consists in 
writing an interface equation for the present model. A 
single-valued one dimensional interface y — h{x, t) is as- 
sumed, and the Green's function of the problem is defined 
through the relation 



(8) 



for the range ~oo < x,x' < oo, < y,y' < oo, with 
Dirichlet boundary conditions. The half-plane must be 
used, since the presence of an "infinite reservoir" at posi- 
tion y = breaks the translational symmetry in y. With- 
out any loss of generalit y, an = is set from now on. 
The standard procedure |2|,^,Q, exposed in Appendix 
may then be followed to obtain the integro-differential 
equation of motion 



dx' G{x,h{x,t)\x',h{x',t)) 



dh{x\t) 
di 



r]{x, h(x, t)) + ct/C 
with the half-plane Green's function 

G{x, y\x ,y) = — \n 

47r (x — x'Y + (y + y ) 



(9) 



(10) 



The quenched noise is written as ri{x,h) = J dy(pQ{y — 
h(x, t)) a{x, y) ~ 2a{x, h) in the sharp interface limit. The 
Gibbs-Thomson boundary condition, fj,\int ~ ^ can be im- 
mediately obtained from Eq. (||) in the limit h = since 
r] is the chemical potential at the interface. Analogous 
non-local equations will arise in the context of directional 
solidification, pattern selection in Laplacian fluid flow [|4| 
and step growth . The novel features here are the bro- 
ken translational invariance and the presence of quenched 
noise. The interface fluctuations are thus intimately cou- 
pled to both the average position and the average velocity 
of the interface, a result that comes out self-consistently 
from the model. This is quite different from local types 
of equations or models. It should be particularly noted 
that the presence of a conservation law does not result 
in a "conserved" interfacial equation. Likewise, nonlinear 
equations with long range kernels |^ do not apply to the 
situation encountered here. 

The difference between local models and Eq. (^ be- 
comes even clearer if the interface is linearised in small 
deviations around the mean interface position H(t) to ob- 
tain 



(l-e-2|^l«)+|fc| 



-2\k\H 



+ \k\Hhk [l + e- 

\k\{Mt)}k~'jk^hk) , (11) 

where hk are the Fourier components of h and H — Hq 
is the average interface position. Note that the interface 
configuration enters the disorder term in Eq. (^Tj) in a 
fundamentally nonlinear way. 



{V{t)}k 



-ikx 



r]{x,h{x,t)). 



(12) 



This equation already yields important information for 
the roughening process. For example, the limit /c — > re- 
produces the slowing down of the front as given by Eq. (^. 
It also reveals the different length scales in the problem. 
The average height of the interface separates two regimes 
of lateral scale. For kH ^ 1, Eq. (O) becomes 



dhk 
~dt 



1 dH 
H~dt' 



G 

2H 



hk + TT^k^hk 



2H 



(13) 



and in the opposite short scale limit, kH ^ 1 



^ + |fc| '^hu + a\k\^K = k{,mk- (14) 
dt dt 

The average interface height H is thus also a lateral 
length scale. If two points are separated by a distance r ^ 
H , they are not connected through the bulk of the sys- 
tem and receive liquid from the reservoir independently. 
In that sense, the dynamics of the interface on larger scales 
is local, although the slowing down of the interface is inher- 
ently a non-local phenomenon, refiecting liquid transport 
through the medium. In the opposite limit oi r ^ H, the 
two interface points will be coupled through the bulk and 
compete for liquid coming from the same region behind 
the front. The dynamics of the interface then becomes 
fully non-local. 
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However, in both limits the damping terms induce a 
separating length scale = {'^/HY^'^ = (cri?/a)^/^. For 
3> 1 the fluctuations of the interface are damped due 
to the hue tension cr, while for ^xfc ^ 1, it is due to flow 
from the reservoir. By this mechanism, it is expected that 
the front is smoother on length scales larger than ^x as 
compared to smaller scales. 

This length scale is closely related to the MuUins-Se- 
kerka instability of driven Laplacian fronts js^, although 
the situation is reversed here. Because fluid is transported 
towards the front from behind, advanced (retarded) parts 
of the interface receive less (more) mass than the aver- 
age and the front is stabilised at long length scales. This 
result can intuitively be understood as follows. Due to 
the Gibbs-Thomson effect a local "bulge" of vertical ex- 
tent W and lateral size ^ alters the chemical potential by 
Aji ~ aW/^"^. On the other hand, the average gradient 
in PL in the bulk liquid induces a difference ~ aW/H 
across a vertical distance W . These two differences bal- 
ance each other at a length given by 

Cx ^ ^<jH/a. (15) 

The length scale ^x is a static quantity, but in a rising 
liquid column it becomes dynamical through the time de- 
pendence of H{t), i.e., Cx = ?x(i) [H{t)Y/^ ~ 
must also be noted that it is not a truly dynamical cor- 
relation length in the sense of kinetic roughening theories 

However, it is a time dependent upper cutoff for corre- 
lated fluctuations increasing with time, and therefore can 
be interpreted as a dynamic correlation length. 

To draw any further analytical conclusions from Eq. 
( pi] ) is extremely problematic due to the difficulties en- 
countered with the quenched noise {77 1^ (which are anal- 
ysed e.g. in [Q). Furthermore, Eq. ( prl ) is a linear ap- 
proximation to Eq. (H). Although the length scales come 
out correctly, it cannot be expected a priori that the cor- 
rect scaling properties of the interface will be obtained. 
An analysis along the lines of Ref. may prove insufff- 
cient, and simplifled treatments, such as those attempted 
in Ref. are inadequate. 

3 Numerical Analysis 

The interface fluctuations in the presence of quenched 
disorder were analysed by numerical integrations of Eq. 
(^ and Eq. (||). The position of the interface h{x^t) at 
each X was deflned by the zero of the phase field, i.e. by 
(t>{x, h{x,t)) = determined by linear interpolation be- 
tween the points of the numerical grid. Overhangs, ap- 
pearing for strong disorder but otherwise absent were ig- 
nored by taking the lowest or highest zero of (j) above a 
given foot point x. No quantitative differences were seen 
between these two choices. 

The disorder a(x) is an independently distributed ran- 
dom variable on each grid point, with mean a, standard 
deviation Aa and chosen from different types of distri- 
butions (gaussian, uniform on a finite interval, and expo- 
nential). Without loss of generality, the lower boundary 



condition is chosen such that p{x,y ~ 0) = 0, leading to 
(j>{x,y = 0) = (fiQ, with (j)Q the solution of —(j>o + (j>Q = a. 

To evaluate any scaling behaviour the first quantities 
of interest are the total width of the front, 

W\t) = {{h{x,t)-H{t)r), (16) 

and the related spatial two-point correlation functions of 
the qth moments 

Gq{r,t)^(\h{x + r,t)~h{x,t)\iy/'^. (17) 

The case g — 2 is directly related to the structure factor 
S{k,t) = {hk{t)h^f^{t)) . In the above equations the brack- 
ets denote an average over different realisations of a, and 
the overbar a spatial average over the system. In presence 
of a stationary state also temporal averages can be taken. 

The standard Family- Vicsek scaling assumption rests 
on a dynamical correlation length ~ t^^^ where z is the 
dynamical exponent, related to the decay of fluctuations 
along the interface. The maximal value it can attain is the 
system size L at which point the interface is said to be 
in a "saturated" stationary state. The two-point correla- 
tion function then has a scaling form G2{r, t) = r^f{r/^t) 
where f{u) = const, for u <^ 1 and f{u) ^ for 
w 3> 1, a form which introduces the roughness exponent 
X and deflnes the associated growth exponent (3 = xl 
The structure factor has a corresponding scaling form, 
S{k,t) = s{k^t)/k^~^'^^ with the scaling function s{u) con- 
stant for u 1 and s{u) ~ for u ^ 1. 

This picture may turn out to be incomplete or even 
wrong for the following reasons. First, the structure fac- 
tor may contain an explicit time dependence besides ^t, 
S{k,t) ~ t^'^ s{k£_t) / k^^"^^ , which is sometimes referred to 
as intrinsic anomalous scaling [^ . Second, if the inter- 
face is superrough, a case characterised by x > 1, then 
G2{r,t) ^ Q {f / £,1)^^°" with a local exponent xioc = 1, 
since, by construction, G2{r) cannot increase faster than 
r [Q. In contrast to the standard Family- Vicsek picture, 
in both these cases the correlation function for r < 
do not saturate as long as ^4 increases [0. This can be 
parametrised by the scaling of the local slopes G2{r = 
l,t) ~ t(x-xioc)/z^ Third, the lateral £,t and the vertical 
scale W may not be enough to characterise the interface 
fluctuations, and different moments of Gq may possess dif- 
ferent scaling exponents Gq{r) ^ (see e.g. pl|). 

The interface scaling behaviour may also be observed 
in the temporal correlation functions 

Gq{t) = {\h{x, t + s)~H{t + s)- h{x, s) + i7(s)|9)i/«, 

(18) 

which increase as Cq{t) ^ at short time differences 
t. Of course, this definition makes only sense in a steady 
state, under time-translational invariance. It therefore ap- 
plies to the analysis of Eq. (|^), where the average interface 
height is kept fixed by pulling down the paper at constant 
velocity v. 
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3.1 Freely rising fronts 

In this Subsection, the classic imbibition experiment is 
considered. A hquid front is allowed to invade the porous 
medium starting from a reservoir placed at y — 0. Succes- 
sive configurations obtained from numerical integration of 
Eq. are presented in Fig. ||. The time difference be- 
tween the curves is constant {At = 10"^), and the slowing 
down of the interface positions becomes apparent from the 
fact that they lie closer together the higher the front gets. 

In the presence of quenched disorder a{x) the total 
width of the interface increases as a power of time. Fig. || 
shows that W{t) ~ with f3 ~ 0.3. In the same graph the 
progression of the average interface height H is seen to fol- 
low Washburn's behaviour as expected from the analytic 
arguments. 

It turns out to be impossible to determine the global 
roughness exponent x from its definition in terms of the 
dependence of the saturated width on the lateral system 
size. Instead, the structure factor S{k, t) is used. It is plot- 
ted on Fig. H for a system of size L = 256, with a = 0.2 
and Aa = 0.2 at various times 10'^ < t < 5x 10"* (in the 
dimensionless units of Eq. (||)), corresponding to heights 
ranging from 20 < H < 100. Although it is difficult to 
obtain good statistics for this quantity, immediately ap- 
parent is a strong power law decay, S{k) ~ l/k^~^^^, with 
a global roughness exponent x ~ 1.25, and a crossover 
to a plateau corresponding to distances larger than the 
time-dependent correlation length ^x(^)- 

The form of the structure factor indicates the presence 
of anomalous scaling; with x > 1 the surface is superrough 
||l[ and the spatial correlation functions will reflect this 
fact. The level of S{k, t) in the region of the power decay 
also seems to increase with time, i.e. S{k,t) ^ t^/k^'^^^, 
with 9 « 0.05, which could indicate the presence of intrin- 
sic anomalous scaling A clear identification of this 
regime is however quite difficult, due to the very slow in- 
crease in time and to poor statistics. 

The anomalous form of the scaling is most visible in the 
two-point spatial correlation function G2 {r, t) as shown 
on Fig. H for the same data as for the structure factor. 
The correlated roughness of the interface is visible up to 
a length scale rmax ^ Cx(i), and the average mean step 
height G{r = l,i) - Cr'"" ~ t(x-x.oc)/4 |||9j^ 

The average driving force a affects the scaling of the 
structure factor and correlation length only through the 
correlation length {t) ~ arid the total strength 

of the noise Aa only influences the amplitude of the pref- 
actor of the correlation function. The correlation function 
may be fitted to the function 



sequence of anomalous scaling q Likewise, the structure 
factor may be described by the scaling form 

Sik,t)^'^§^, (20) 

where the scaling function s(x) is constant for x ^ 1 
and s{x) ~ x^^"^^ for a; <C 1. The scaling behaviour can 
be seen in the inset of Fig. ||, where the scaled form of 
the correlation function is shown for a single system of 
L — 256, with a = Aa = 0.2 at various times and in Fig. 
^ where the scaled correlation function is now shown for 
systems of similar lateral extent but different values of the 
driving force and strength of disorder. Within this scaling 
picture, the early time development of the width follows 
W{t) ~ ~ ^'^^^ = yielding a growth exponent (3 = 
x/4 ft! 0.31 in good agreement with the direct numerical 
estimate. 

3.2 Fronts at Fixed Height 

The results presented in the last section indicate a scaling 
picture to be valid in the freely rising case. It is however 
difficult to obtain sufficient statistics and larger samples 
of quenched disorder replicae are necessary to get accurate 
data. This difficulty can be overcome by considering the 
stationary interface, as described by Eq. (^). In this case, 
the interface fluctuations reach a steady state and the var- 
ious correlation functions can be obtained with greater ac- 
curacy. This setup also allows the investigation of height 
difference temporal correlation functions and permits a 
comparison with the experimental results of Horvath and 
Stanley |l|. 

Equation (j^) was integrated numerically for different 
values of mean height H = a/2v (see Eq. (|)). Differ- 
ent lateral system sizes were used between L — 12 and 
i = 400. The total vertical extent of the lattice was taken 
to be about 50 length units higher than the interface, to 
prevent any influence of the upper boundary. The mean 
value a — 0.2 in all cases, but different distributions were 
used: (i) a uniform distribution with mean a = 0.2, on 
the range [0, 0.4) (standard deviation Aa = 0.07), (ii) the 
same, but on the range [0.1,0.3), {Aa = 0.03 ), (iii) an 
exponential distribution with average 0.2 {Aa — 0.2). Dif- 
ferent numbers of configurations were used in taking the 
averages, from 10 in the largest systems {L = 2H = 400) 
to 100 in the smallest (L = 2i/ = 50). Even in the largest 
systems, saturation of the interfacial fluctuations became 
apparent after times t ~ 2 x lO'*. All systems were inte- 
grated up to i = 10^, with the interface configurations 
extracted at time intervals At = 100. 



G2{r,t)^Aae^g{r/S,^), 



(19) 3.2.1 Spatial correlations 



with X = 1-25 and a scaling function g{x) 



with f{x) ~ x~-'<^'°'= for a: ^ 1, and approaching a constant 
for a; <C 1. There seems to be no simple explanation why 
Aa enters Eq. ( p^ ) in a linear way even beyond the linear 
approximation of Eqs. (|l^) and (14) where it is easy to 
see. The local scaling exponent xioc — 0.9 is a direct con- 



Along with the height H , the length scale ^x also remains 
fixed in this setup, since it is related to the driving ve- 
locity through ^x 



V 



-1/2 H^/"^. Thus, contrary to the 

^ In principle xioc ~ 1 for a superrough interface. However, 
here the scaling behavior in ^(A:, t) appears only over a rela- 
tively short range and finite size effects are pronounced. 
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standard picture of kinetic roughening, the saturation of 
the interface is not necessarily determined by the total 
lateral extent L of the system. Here, the correlations sat- 
urate at either the system size L or the correlation length 
J whichever is smaller. Fig. |^ shows data obtained for 
= 50 on system sizes L = 12, 25, 50, 100, 200 and 400. 
The curves for L = 100, 200 and 400 collapse indicating 
that the fluctuations are bounded by the L independent 
length scale (50<^x < 100 for this particular case). 

The structure factor S{k,H), shown in Fig. ^, also 
has a pronounced power law decay A:"^^-^"'"^-' with a global 
roughness exponent %= 1.25. Again, the interface is su- 
perrough. As in the freely rising case, there also seems to 
be a very weak intrinsic anomaly in the sense of |^9| , i.e. 
the prefactor of S{k, H) in the power law region depends 
on H. The data are consistent with for O<0<O.l, but 
not accurate enough to draw any firm conclusion here. 

Provided that < L the spatial correlation function 
G2 {r, H) follows the scaling form 

G2{r,H) = Aav-"^'^ g{rv^'^). (21) 

with the scaling function g{x) defined as in Eq. (|l|). For 
aU different setups, L = 50 to 400, = 25 to 200, and 
all three choices for the disorder a(x), G2{r,H) is shown 
rescaled according to Eq. (pT| ) in Fig. For small dis- 
tances r < ^x the spatial correlation function is of the 
form G2{r) ~ r^x yielding a local roughness exponent 
Yioc~ 1 and a height difference at fixed r growing as ^x~^ 
§ (see Fig. I). 

Finally, the different moments g = 2,4 and 6 of the 
correlation functions Gq{r) can be compared, as shown in 
the inset of Fig. ^. All moments have a local exponent 
Xioc ~ 0.95. The global exponent Xq can in principle be 
obtained from the short distance scaling oi Gq{r ~ 1) ^ 

S^" ^ . The present data point to a similar value Xq = 1-25 
for all q's but are however too noisy to draw any definite 
conclusion. Thus, the interface is probably truly self-affine 
up to the crossover scale — be it the system size L or the 
saturation length ^x- 

3.2.2 Temporal correlations 

Before the interface fluctuations reach the steady state 
they are governed by an increasing dynamical correlation 
length. Starting from a flat front h{x,t = 0) = it is 
observed to grow roughly as t^^^ and then approach 
^x ■ However, this behaviour could not be analysed in much 
detail, due to the short time range of the initial power 
law and to poorer statistics (averaging over time is not 
possible). The insets in Fig. shows S{k,t) approaching 
the saturated regime for the large system, L — 2H — 400. 

Next, the correlation functions G2{t), shown in Fig. ^ 
for different heights are compared, as in the experiments of 
Horvath and Stanley |jl^. In both cases, the crossover time 
ts between the power law regime of C2{t) and saturation 
increases with H. This is true for the level of saturation (or 
the width W) as well. At early times i <C is , the absolute 
value of C2(i) decreases with H. 



The data can be related to a scaling function of the 
form 

G2{t)^m/'f{t/H^/'), (22) 

with the scaling function f{x) ~ x^^ for x <^ 1 and const. 
for a; S> 1 with a (genuine) dynamical exponent 2; « 2 and 
the effective slope 02 ~ 0.85. Although the exact value 
of these exponents is difficult to establish, this form is 
however valid for all L provided that ^x < L. 

The different moments Gq{t) for g = 2,4 and 6 are 
shown in Fig. |l^. They clearly have different behaviour. 
The early time logarithmic slopes of the higher moments 
(Pq) decrease with q, as shown in the inset of Fig. ^ For 
the higher moments, the effective exponents P4 ^ 0.76 
and Pe ~ 0.69. Such multiscaling has been observed in 
cases connected with the existence of avalanches in the 
interface dynamics It is clear from Fig. ||that similar 
avalanche type of motion exists here as well, but only up 
to the vertical length scale W and lateral size ^x ■ Because 
of this reason the quantitative characterization of such 
avalanches is beyond the scope of the present work. 



4 Discussion 

4.1 Temporal Scaling of the Interface and Relation to 
Experiments 

It is interesting to first compare the results for the station- 
ary fronts with those of the freely rising fronts. Both cases 
are governed by the same height dependent length scale 
^x- In Fig. ^ the spatial correlation function G2{r, H) 
(fixed H at saturation for the stationary case) and G2(r, tn) 
(freely rising case at times tn = H'^/a when the average 
height has reached H) are shown for various values of H. 
There is a complete equivalence between the interfacial 
fiuctuations at an instantaneous height H{t) and the sat- 
urated fluctuations of a stationary interface. In both cases, 
the range of correlated roughness is determined by the 
same value ^x ■ The length scale ^x is thus conceptually 
different from the intrinsic time dependent lateral correla- 
tion length commonly found in models of kinetic rough- 
ening. Here ^x merely fixes the maximum range of corre- 
lated roughness. Such "quasi-stationarity" of the moving 
front can only occur provided that the "natural" dynam- 
ical exponent z < 4, so that the interface fluctuations can 
always catch up instantly with the available area of cor- 
relation. 

The model can also help to interpret the experiments 
performed by Horvath and Stanley on the stationary inter- 
face ll^, since it yields qualitatively similar results: The 
exponent /32 is constant for all driving forces (and also f3q 
for the higher moments considered) , the level of saturation 
of C2{t) increases as H increases, while the amplitude of 
the early time power law behaviour decreases. 

Still, there is a quantitative difference between the ex- 
periments of ]l3| and the numerical results. In the present 
work, P2 =0.85 in contrast to the experimental value 0.56. 
Clearly one would not expect coincidence, since already 
the average front velocity behaves differently: The model 
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shows "pure" Washburn behaviour, dH/dt^l/H, whereas 
in (13) dH/dt ~ H''^-^. This deviation from Washburn's 
law occurs in many paper imbibition experiments [ p4| , and 
there is no reason why its origin should not affect the in- 
trinsic fluctuation dynamics as well, leading to different 
values of Pq in the time correlation functions. 

There is however a clear contradiction of the results 
with the scaling form for the time-correlation function sug- 
gested in [O, 



(23) 



with the relation V = dH/dt ^ H and the scaling 
function C{u) such that C2{t) ^ f^V^* if u < 1 and 
C2{t) ~ L^V~^^ in the opposite limit. Equation ( p3| ) does 
not include an explicit lateral length scale ^x, although 
some power of y^'+^i- may obviously play such a role. This 
scaling form is however in contradiction with the present 
results in two ways. First, provided that > L, the total 
lateral system size does not play any role in the scaling 
of the interface. Secondly, no common scaling form inter- 
polating from large systems (fluctuations up to scale ^x ) 
to small systems (dominated by L) has been found. The 
saturation of C2{t) is sharper for small L, which is visible 
in Fig. |ll|. 

In the experiments only one lateral system size was 
used, so the role of the system size L could not be assessed. 
Unfortunately, no information on the spatial scaling of 
fluctuations was presented so that no comparison can be 
done to the present work. Likewise, the higher moments of 
the temporal correlation function were not measured and 
the presence of avalanches could not be inferred. 



4.2 Avalanches, Pinning and Roughness of the 
Interface 

As one can see from Figure ^ the interface motion indicates 
the presence of "avalanches" as in usual models describing 
depinning transitions j^. Here the behavior is however 
somewhat different since the flow of liquid in the average 
tends towards regions of lower chemical potential. In some 
way, this is analogous to 'self-organised' interface models 
in which the interface is driven at the point where the 
force is the largest HQQ. 

The question is now whether it is possible to under- 
stand the observed exponents — Xj Pg 6tc. — in terms of 
an avalanche description. Were this to be true, the pres- 
ence of the conservation law would only be felt through 
the correlation length ^x , which would limit the avalanche 
area to ^x''^ with an exponential cut-off on the avalanche 
size distribution (see for a discussion of a related ex- 
periment). It turns out that the multiscaling observed 
here (Pq) differs in a crucial way from that obtained for 
avalanches in self-organised depinning. Leschhorn and Tang 
obtain an almost trivial multifractal spectrum for the 
Pq^s assuming that a local dynamical exponent zioc and 
a global roughness exponent can be defined. Dynamical 
scaling in this sense is however absent in the present case, 



since the interface dynamics depends directly on the local 
height. 

The observed global roughness exponent x — 1-25 also 
appears in "nonconserved" front propagation through a 
medium with quenched disorder |^,^,^, i.e. for inter- 
faces in the quenched Edwards- Wilkinson (QEW) univer- 
sality class. In that problem a similar value of x is observed 
close to criticality, whereas a cross-over to thermal EW 
(the massless Gaussian field, with x — 0.5) takes place in 
the moving phase. The imbibition model may then present 
a similar behaviour due to the continuous slowing down 
of the interface, with parts of the interface pinned (thus 
approaching the depinning transition from above). As in 
the QEW class, the value x = 1-25 means that the width 
increases faster than the correlation length. The "local 
slopes", G2{r = l,t) ~ ^x~'^'°''(^) diverge with time and 
the local roughness exponent xioc ~ 0.9 . . . 0.95. For QEW 
models this is true on any scale — up to saturation in a 
finite system — but in our case the behavior can hold only 
up to ~ ^x , because the removal of overhangs occur 
naturally in a phase field formulation. If W > ^ X over- 
hangs in neighboring bulges of the interface can merge, 
always keeping ^x of the order of W. 

A possible mechanism by which the roughness might 
change is if the interface becomes completely pinned. In 
Appendix it is shown that this may be achieved through 
the inclusion of either evaporation or gravity. In such a 
case, new length scales come into play, and the roughness 
exponent is most likely changed. 

Although the pinning effects of evaporation have not 
been explicitly considered in this work, preliminary calcu- 
lations show that evaporation, characterised by an evap- 
oration rate e introduced in the Appendix pins the in- 
terface at a height Hp ~ (a/e)-'^/^. At this pinning height, 
a new correlation length emerges. For weak evapora- 
tion, defined as e <C c? ja"^ , the length ~ (cr/e)^/"^ S> 
^x{Hp). For strong evaporation, e ^ a^/a'^, this length 
is ^ i'^/^HpY^^ ^ (,x{Hp). This question, as well as 
the appearance of ^x itself, has not been addressed by 
previous models and experiments [ p^ . 

Within this region, there should be a crossover from 
the superrough interface with x = 1-25 to a pinning regime. 
The roughness exponent on scales below the correlation 
length should then be determined solely by the local dis- 
order configuration, i.e. be related to e.g. directed perco- 
lation depinning. 



4.3 Hydrodynamical description of Imbibition 

The model presented in Section |^ is purely diffusive and 
does not include any hydrodynamical modes. In princi- 
ple, these could be incorporated, in a coarse-grained sense, 
along the lines of Refs. Q and Q by coupling the phase 
field to a velocity field described by the Navier-Stokes 
equations. Unfortunately, many problems need to be re- 
solved before such an approach is taken, the main one 
consisting in establishing the role of hydrodynamics itself 
[ p4| . Contrary to other porous media like fractured rock 
or Hele-Shaw cell filled with glass beads, the paper matrix 
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used in most imbibition experiments is not inert but of- 
ten interacts strongly with the invading fluid through fiber 
svifelhng. Another comphcation arises from the transport 
of fluid through the paper. It is not at all obvious that 
an homogeneous pressure can be defined throughout the 
volume occupied by the fluid. In that sense, the model 
defined by Eq. (H) and Eq. is the minimal model that 
includes a conservation law and, in conjunction with the 
appropriate set of boundary conditions, reproduces the 
experimental characteristic of imbibition. 

Even though hydrodynamics is absent there is how- 
ever of course a strong similarity between the model and 
the standard description of flow in porous media based 
on Darcy's Law for an incompressible fluid. In this de- 
scription, the normal velocity of the interface is related 
to the gradient in the pressure field P(x) by the per- 
meability K as Vn — —K dnP{x,y = h{x,t)). The pres- 
sure is determined from Laplace's equation, together with 
P{x,y = 0) = Pq, the atmospheric pressure, and a Gibbs- 
Thomson boundary condition at the interface: P{x, h) = 
Po—Pc{x, h{x, t))—jV^hJx, t), which introduces a coarsed- 
grained surface tension and capillary pressure Pc aris- 
ing from the microscopic menisci at the fluid-gas interface. 
Working to linear order in the small fluctuations of the 
interface, it is straightforward to find the pressure field, 
defined for y < h{x,t), 

P{x, y) = P„ - P,j^ e"^-^ sinh(fcy) P^ , (24) 

k 

with the coefficient 

^^^-iml^(G^^-|)'^^-^^(^^«'^0' 

(25) 

where Pc represents the average capillary pressure and 
Pc{k 7^ 0, H) the fluctuations around it. It is then a sim- 
ple matter to derive the interface equations, Eq. ( [T^ ) and 
Eq. ( p^ in the appropriate limits. The derivation of the 
length scale may be transposed directly to the fluid 
imbibition case. This is essentially equivalent to the ap- 
proach of Ref. (see also [0), although spontaneous 
imbibition requires a special treatment of the boundary 
conditions, absent in their work. 

Our model is based on a constant mobility, which in the 
general case should be replaced by dtcj) = VA/(x)V^(x, t). 
In spite of its simpleness, it is quite reasonable to ask 
whether the quenched randomness should not be included 
in the mobility, to model an effective quenched perme- 
ability. In most cases of forced fluid flow in bulk random 
media, this is where the non-uniformities are most rele- 
vant iQ. One point must however be emphasised. In any 
imbibition experiment designed for this purpose, the flow 
will never be large, and we believe that the random capil- 
lary forces will have the dominant influence. On the other 
hand, it has recently been shown that the presence of ink, 
or presumably of any other blocking material does make 
a quenched porosity relevant pl| ] , a case which is not con- 
sidered here. 

It should also be pointed out that the field (p does not 
represent a real fluid density. In particular, any "air bub- 



bles" (i.e. connected regions of value (p = —1) trapped 
behind the front will eventually dissolve. This is of course 
highly unrealistic for bulk porous media, but may be ap- 
propriate for thin porous media, where air can escape 
through surface pores. 

5 Conclusion 

In conclusion, a simple phase field model for the invasion 
of a liquid into a disordered system has been introduced. 
Liquid conservation is explicitly included. Some basic fea- 
tures observed in imbibition experiments are reproduced. 
Of course, it cannot account for many phenomena of paper 
wetting or invasion into porous media, which are briefly 
discussed in Section ^ and treated in more detail in the 
previous paper Q . 

In numerical simulations, a superrough interface is found, 
with anomalous scaling due to a global roughness expo- 
nent X — 1-25 > 1. A hardly discernible inherent anomaly 
of the structure factor may also be present. The extent 
of the spatial fluctuations of the interface are bounded by 
a length scale ~ i^'H/cxY^'^ , both in the freely rising 
and stationary imbibition setup. Interface fluctuations of 
a rising front are quasistationary, in the sense that at any 
time t they are the same as in a stationary front kept at 
average height H = H{t). The temporal fluctuations show 
multiscaling, which indicates motion by avalanches. The 
length and time scales where fluctuations saturate can be 
understood by simple dimensional considerations. 

In relation to experiments, according to this analysis it 
is highly desirable to have fluids with small capillary pres- 
sure and high surface tension in order to obtain scaling 
over a large spatial regime. It would also seem appropri- 
ate to use organic liquids in experiments done with paper. 
These have minimal chemical interaction with the con- 
stituent fibers, and a simplified Washburn description of 
imbibition may be applicable. Another option is to use 
deionised water again with the goal of reducing the 
fiber-liquid interaction. Above all, the main conclusion of 
the present work is that the macroscopic behaviour of the 
average interface position (i.e. H(t) in a freely rising study, 
or H{v) if the interface is stationary) is crucial to an un- 
derstanding of the microscopic fluctuations of the inter- 
face, since it controls the range over which scaling can be 
observed. 
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A Refining tlie Phase Field Model 



and the equation of motion for (!>tot is 



In this appendix, a generalised imbibition model is intro- 
duced. It rests on the following dynamical equation for the 
dimensionfuU phase field 4>{y, t) 



dT dy 54) 

with the free energy functional 

dv 



1 



^=2 



a(j) 



(r,r)) 
(26) 

(27) 



and 4>e — (r/w)^/^. For (5 = and e = 0, as well as for con- 
stant mobility M{(f)) — M, this reduces to the model de- 
scribed in the introduction. When e ^ 0, a non-conserving 
term is introduced in the equation of motion, which, in a 
first approximation describes evaporation of liquid (the 
(j) = (j)e phase), at a rate 2(j)2e and proportional to the 
total area covered by the fluid. The convective term is 
included to describe gravity. Although gravity can only 
be introduced properly through an hydrodynamical field 
1^,^, it is shown below that this term reproduces the 
correct equation of motion for the average position of the 
imbibition front. 

Assuming a constant mobility M, these equations can 
be put in a dimensionless form by defining 



r/C; 



a = 
G 



I 1/2 ~ 

ll/2 _G . 
J M ' 



(28) 



_ 



My2 



The ratio ( = (ft/r)^/^ determines the width of the inter- 
faces between different phases. 

For now, let us concentrate on the case G = e = 0. In 
one dimension, with a being constant, the chemical poten- 
tial obeys a Laplace equation in both wet and dry phases 
(actually there are minor corrections) with boundary con- 
ditions, fi{y = 0) = and /Lt(y ~ H) = /i(y = Ly) — —a, 
where H is the position of the wet /dry interface and Ly 
is the length of the paper. Thus 



which implies that 



-ay IB, 
-a. 



0o + (1 - <i>o) y/H, 



-1, 



ify <H: 
iiy>H- 



iiy>H- 



(29) 



(30) 



where (/jq > 1 is the solution oi —4> + cf) = a. In a first 
approximation, both (p and jj, are linear functions of y for 
y < H. The total amount of concentration is then 



4{y,t)dy^-[<i>o + m{t) 



(31) 



d<^t 



idH 



dt 



dt 



dy 



^dy, 



dH{t) 
dt 



3 + 00 \H{t)) 



(32) 



(33) 



an equation similar to Washburn's result. Since 0o ~ 1 + 
0(a), this is actually Eq. (|^) of Section |^. In presence of 
gravity and/or evaporation, the solution is more involved. 
As a first approximation, the difference between (f)Q and 
unity is neglected and the Poisson equation for the chem- 
ical potential is considered, 



dy^ 



e = 0. 



for y < H(t) and boundary conditions fj,{y 
and fi{y = H{t)) ^ —a. Again, fi{y > H{t)) = 
solution of Eq. (U) is 



ay 

m 



-ey{y-H{t)) 



(34) 



0) = 
-a. The 



(35) 



Using the same procedure as above, the equation of motion 
of the interface is found to be 



dt 2H{t) 4 ^ ' 



(36) 



The "gravity" term acts exactly as in Washburn's equation 
in presence of gravity, and thus allows to identify G as an 
effective gravity force acting on the interface. As far as 
we are aware, no detailed studies of fluid propagation in 
a thin porous medium with evaporation has been done. 

Non-zero values of G or e will eventually stop the in- 
terface at an equilibrium height Heq — a/G ifG 3> e 
or Heq = (a/e)-^/^ if e S> G. There is however a concep- 
tual difference between pinning due to gravity or evap- 
oration. In the former case, the chemical potential is a 
linear function of position, and pinning is determined by 
dnfJ.{Heq) — G whilc, in the latter case, the chemical po- 
tential is quadratic in y and at pinning, dnfJ-iHeq) — 0. 
When both evaporation and pinning are present, the equi- 
librium height is determined by the zero of Eq. ( ]3^ ) . 

In terms of a dimensionfuU interface height H, Eq. (|3 
becomes 



dHjT) 
dT 



1 Ma 

2 (be 



1 



H{t) 



G 



(37) 



where a represents the average value of the disorder. The 
motion of the average interface in the pure Washburn 
case is thus determined by a combination of the mobil- 
ity M and the shift in the average chemical potential 
a. The length scale ~ {uH/aY/'^ where now, a = 
(2\/2/3)(Kr0g)^/^ is the dimensionfuU surface tension. 
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B Projection to an Interface Equation 

To extract the interface equation Eq. in the hmit G = 
and e = 0, the dynamical phase field equation must first 
be inverted with the use of the Green's function defined 
byEq. {§ and Eq. 



Jd^' G(x|x')^^^ = /i(x,i). 



(38) 



It is then convenient to use a local coordinate system (u, s) 
pO{ . The 2-dimensional space is spanned by the the vec- 
tor x(m, s) = X(s) + un(s), where X(s) is a point of the 
interface, fi is a unit vector normal to the interface and s 
is the arc-length coordinate. In terms of the phase field, 
this corresponds to = 0, s) = 0. The time derivative of 
the field then becomes s, t)/dt = Vn{s)d(j)/du where 

Vn{s) is the normal velocity of the interface at position s. 
If the interface (of thickness C = 1 in dimensionless units) 
is much smaller than the typical radii of curvature of the 
interface (the sharp interface limit), the Laplacian term of 
the chemical potential may be expanded such that 



52 92 B 
ou'^ os'^ ou 



(39) 



where /C(s) is the curvature of the interface. For a — Q, 
the one-dimensional kink solution is (j)(u,s) = (/>o(w) = 

tanh(u/V2)- The corrections to this form, represented by 
Eq. ( |30| ) are of order (d/H. To first order thus, ^ 
—a{u, s)—IC(s)d(j)o{u)/du. Still in the sharp interface limit 
^/C <C 1, the derivatives of the kink solutions have prop- 
erties d(j)o{u)/du ^ A(f)6{u) and cr = f du{d(j)o{u) / du)'^ 
where A(j) 2 is the miscibility gap and a is the interface 
tension. Multiplying Eq. (|3^ ) by / du{d(l)Q{u)/du) then 
effectively project the phase field dynamics onto the in- 
terface u = 0. A translation u —^ u + h{s, t) then yields 



ds' G{s, h{s, t)\s', h{s', t))Vn{s') = Tjix, h{x, t)) + aJC. 

Equation (|g|) is then obtained by a further change of co- 
ordinate s ^ X and the relation ds'Vn{s') = dx' dth{x' , t) . 
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Fig. 1. Setup of the imbibition model. The system is defined 
on the plane y > 0, with lateral extent L. The average position 
of the interface is represented by H{t). The chemical potential 
obeys Laplace's Equation in the bulk V^/i = with Gibbs- 
Thomson boundary condition /x = — a(x) — ad^h{x,t) at the 
interface and the imposed value fj,{y = 0) = ao = at the 
bottom. 
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Fig. 2. Front configurations of a rising interface at equal time 
intervals At = 10"^. Their average separation becomes smaller 
as the front slows down due to the conservation law. 




Fig. 3. Increase of the width of the interface as a function 
of time. The data are for a system of lateral extent L — 256, 
with a = 0.2 and Aa = 0.2. The straight dashed line has slope 
P = 0.32. The inset shows the average position of the interface 
H^{t) = at, in agreement with Eq. (g). 




Fig. 4. Structure factor for a system with L = 256, a — 
Aa = 0.2 at different times. The data are for times t = 5 x lO** 
(lowest curve) to t = 5 x 10"* (upper curve) at intervals of 10*. 
The dashed line has a slope —3.5 indicating a global roughness 
exponent x ~ 1.25. A weak intrinsic anomalous scaling, in 
the sense of Ref. ^] may be present, but cannot be clearly 
identified. 



Fig. 5. Spatial correlation functions G2{r,t) for parameters 
identical to those of Fig. ^. The anomalous scaling is indicated 
by the increase of the average mean step height G2{r = l,t) 
and the local exponent Xian,^ 0-9- In the inset, the data are 
rescaled according to Eq. (y_9|). The scaling shows the existence 
of a lateral length scale ^ ~ t^^'^. 

10° I 




10' 10° 10' 

Fig. 6. Plot of the correlation functions G2{r,t) according to 
the scaling form Eq. (^^. The data are taken from systems 
with lateral size L — 256 and gaussian distribution, with pa- 
rameters (i) a = 0.05, Aa = 0.2, at height H = 20, (ii) 
a = 0.1, Aa = 0.2, at height H = 67, (iii) a = 0.2, Aa = 0.2, 
at height H = 77, (iv) a = 0.3, Aa = 0.2, at height H = 102, 
(v) a = 0.2, Aa = 0.1, at height H = 95. 
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Fig. 7. Spatial height difference correlation functions for a 
setup with fixed average height H — 50, and a(x) uniformly 
distributed on the range [0,0.4). The data for L = 100, 200 
and 400 fall together (top curves), while small systems show L 
dependence (L = 50,25 and 12, middle to bottom). 
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Fig. 8. Structure factors S{k) for H = 25, 50, 100, 150, 200 
(bottom to top) and L = 2H, disorder chosen from an ex- 
ponential distribution with a = 0.2. The power law decay is 
proportional to k~^'^, indicating a global roughness exponent 
X = 1.25. The large scale cutoff k* ~ 1/^x decreases with H. 
The inset shows the approach of S{k, t) to the saturated S{k) 
in the system of size L = 2H = 400 and for times t = 2" x 10^ 
with n = 0, 1, ...7. 



Fig. 9. Correlation function G2{r) for the same data as in Fig. 
^. The curves saturate at a length which increases with H, 
together with the "step height" G2{r = 1) . The local scaling 
exponent xioc ~ 0.95, close to the expected value X!oc = l. In 
inset, the higher moments of the correlation function Gq{r,H) 
are shown for a system with L = 200 and H = 50. All moments 
have the same local exponent Xioc,q ~ 1. 
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Fig. 10. Scaling plot of G2{r,v)/{Aa ^^) vs. r/^x with global 
roughness exponent x=l-25 and ^x ~ v~^^^ ~ H^^^. A wide 
range of H , from 25 to 400, and all three forms of disorder are 
used. The data compare well to the scaling relations derived in 
Eq. (§). 
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Fig. 11. Correlation functions C2{t) for L = 400 and H — Fig. 13. Comparison of the correlation function G'zij',!!) in 

50, 100, 200 (solid curves). The crossover time ts and the level the stationary setup (solid lines) at heights H = 25, 50 and 100 

of saturation increase with H. At short times Czit) decreases to G2{r,tH) in the freely rising case (dashed lines) at corre- 

with H. In inset, C2{t) for L = 50 and _ff = 100 (dashed line) is spending times tn = H^/a. There is a complete equivalence 

compared to C2{t) for L = 400 and H=100 (solid curve). The between both situations, 
transition to the saturated regime is sharper for smaller L. 
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Fig. 12. Correlation functions Cq{t) for g = 2,4 and 6 for 

systems of size L = 200 and H = 50, such that L > ^x- Each 
moment increases with a different exponent Pq. 



